%function [rects] = applyLinearSVMToImage(im, W, SVMthresh, posRECTexamples)
% im = image
% W = linear SVM weight vector
% SVM threshold (try 0.5?)
% posRECTexamples: examples of positive rectangles, used to find reasonable
% sizes.
function [rects rr] = applyLinearSVMToImage(im, W, SVMthresh, posRECTexamples)

% idea, need to scale image so that *every* original rectangle size
% creates 
R = posRECTexamples;

% make exponentially distributed set of image aspect ratios to try
SA = exp(linspace(log(min(R(:,3))),log(max(R(:,3))), 6));
SB = exp(linspace(log(min(R(:,4))),log(max(R(:,4))), 6));
aspectMin = min(R(:,3)./R(:,4));
aspectMax = max(R(:,3)./R(:,4));

[mA mB] = meshgrid(SA,SB);
mA = mA(:);
mB = mB(:);

%parfor mx = 1:length(mA)   % i've flatten the mA, mB mesh so there is only one (longer) loop
fprintf('computing hog');
RSALL = zeros([size(im,1) size(im,2) length(mA)]);
for mx = 1:length(mA)   % i've flatten the mA, mB mesh so there is only one (longer) loop  
    sA = mA(mx);
    sB = mB(mx);
    if sA / sB > aspectMin & sA/sB < aspectMax
        im2 = imresize(im, [size(im,1).*100/sB, size(im,2).*50/sA]);
        hog = vl_hog(single(im2),10);
        
        rectScores = zeros(size(hog,1), size(hog,2));
        for ix = 1:size(hog,1) - 10
            for jx = 1:size(hog,2) - 5
                hogRECT = hog(ix:ix+9,jx:jx+4,:);
                hogRECT = hogRECT(:);
                hogRECT(end+1) = 1;
                %  hogRECT = hogRECT ./ sqrt(sum(hogRECT.^2));
                rectScores(ix,jx) = W' * hogRECT;
            end
        end
        rSnorm = imresize(rectScores,size(im(:,:,1)));
        RSALL(:,:,mx) = rSnorm;
%          imagesc(rectScores);
%          title([sA sB]);
%          drawnow;
        fprintf('.');
        aspectOK(mx) = 1;
    else
        aspectOK(mx) = 0;
    end
end
fprintf('\n');

% find best score at each pixel location
RSALL(:,:,aspectOK==0) = [];
mA(aspectOK==0) = [];
mB(aspectOK==0) = [];
[rsMax rsIdx] = max(RSALL,[],3);

imMax = ordfilt2(rsMax,169,ones(13,13));
rectScores = imMax;
rectScores(rsMax<imMax) = 0;

% scaling to make all dvs [0,1]
rectScores = rectScores - repmat(min(min(rectScores)), size(rectScores));
rectScores = rectScores./ max(max(rectScores));

% crappy threshold
rectScores(rsMax < SVMthresh) = 0;

% toppercent = round(0.1*length(rectScores));
% % crappy threshold version 2:
% rr = rectScores(rectScores>0);
% if length(rr)>toppercent
%     rr = sort(rr);
%     newThresh = rr(end-toppercent);
%     rectScores(rsMax<newThresh) = 0;
% end


% with ties, could be multiple points next to each other.
% label the connected components and just get back the centroids as a cheap
% fix.
BW = bwlabel(rectScores);
regions  = regionprops(BW, rectScores, 'centroid','MeanIntensity');

% get all the scores for the top 21 regions
rr = cat(1,regions.MeanIntensity);
[rval ridx] = sort(rr);
% re-order them by score...
regions = regions(ridx);

topLefts = cat(1,regions(:).Centroid);
topLefts = round(topLefts); 

if isempty(topLefts)
    rects = [];
else
    rects = topLefts(:,[1 2]);
    for rx = 1:size(rects,1);
        sizeIDX = rsIdx(rects(rx,2), rects(rx,1));
        rects(rx,3) = mA(sizeIDX);
        rects(rx,4) = mB(sizeIDX);
    end
end

% show em, to debug?
% for ix = 1:size(rects,1), rectangle('Position',rects(ix,:)); end

